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Abstract. In this paper we describe a quantum algorithm to solve sparse sys- 
tems of nonlinear differential equations whose nonlinear terms are polynomials. 
The algorithm is nondeterministic and its expected resource requirements are 
polylogarithmic in the number of variables and exponential in the integration 
time. The best classical algorithm runs in a time scaling linearly with the num- 
ber of variables, so this provides an exponential improvement. The algorithm 
is built on two subroutines: (i) a quantum algorithm to implement a nonlinear 
transformation of the probability amplitudes of an unknown quantum state; 
and (ii) a quantum implementation of Euler's method. 



1. Introduction 

Systems of nonlinear differential equations arise in an astounding number of 
applications across the sciences ranging from engineering, biological systems, and 
mathematics. Their theory occupies a large proportion of the scientific literature 
and many subtle and profound techniques have been developed to solve them. 
The numerical solution of nonlinear ODEs is now a mature and well-established 
topic (see, eg., [17]) and there are many stable and efficient classical algorithms to 
numerically integrate these equations, ranging from the "workhorse" Runge-Kutta 
method, to predictor-corrector methods and implicit methods. 

The discovery of quantum algorithms (see, eg., [13] for an introduction to quan- 
tum computation) has ushered in a new era where previously intractable problems 
can now be theoretically solved efficiently on a quantum computer. Such was the 
optimisation generated by the discovery of the methods such as Shor's factoring 
algorithm that many felt it would be a matter of time before efficient quantum 
algorithms would be found for many natural problems. This optimism has been 
largely diminished by the realisation that efficient quantum algorithms exploit sub- 
tle - and still largely mysterious - properties of delicate quantum superposition. 
Despite these complications several flavours of quantum algorithm have been de- 
veloped exploiting, variously, algebraic structures, symmetries, and geometry (see, 
eg., [5] for a recent review). Arguably some of these algorithms have less practi- 
cal utility than, say, Shor's factoring algorithm, and constitute more of a proof of 
principle. 

Recently an efficient quantum algorithm to solve systems of linear equations 
[8] was described. This promises to allow the solution of, eg., vast engineering 
problems. This result is inspirational in many ways and suggests that quantum 
computers may be good at solving more than linear equations. In this paper we 
investigate this hope in the context of nonlinear ODEs and we present an efficient 
quantum algorithm to integrate large sparse systems of nonlinear ODEs. 
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2. Preliminaries 

In this paper we are concerned with the solutions of a set of n first-order nonlinear 
ODEs whose nonlinear terms are given by n polynomials f a (z) in n variables Zj, 
j = 1, 2, ... , n, over C. That is, we look for solutions z(i) of the simultaneous set 



d Zl (t) 

dt 
dz 2 (t) 

(1) dt 



dz n (t) 

dt 



= fi(z 1 (t),z 2 (t), . . .,z n (t)) 

= f2(zi(t),Z 2 (t), . . -,Z n (t)) 

= f„(zi(t),z 2 (t), . . -,z n (t)), 



subject to the boundary condition z(0) = b. Standard results ensure that a solution 
to this initial value problem exists and is unique (see, eg., [1]). 

We'll mostly only describe the algorithm for quadratic systems (the extension to 
higher degrees is straightforward). Quadratically nonlinear equations can exhibit 
a wide variety of phenomena, including, chaos and anomalous diffusion and clas- 
sical examples include the eponymous Lorenz system and the Orszag-McLaughlin 
dynamical system. 

We encode the variables Zj(t) as the probability amplitudes of a quantum state 
of an (n + 1) -level quantum system: 

n 

(2) |«/>> = -=|0> + -=]>>|j>, 

where, to ensure that the state is normalised, we require Y^j=i \ z j\ 2 = 1- (Our 
constructions ensure that \<j>) always remains normalised.) While it is convenient 
to regard \<j>) as the state of a single (n + l)-level quantum system, when actually 
implementing the algorithm on a quantum computer we'll encode \<j>) as a state of 
log(n) qubits in the natural way. 

We describe our procedure in two stages: (1) a method to effect nonlinear trans- 
formations of the probability amplitudes of \(j>); (2) a quantum algorithm to imple- 
ment Eulcr's method. 



3. A QUANTUM ALGORITHM TO EFFECT A NONLINEAR TRANSFORMATION OF 

THE AMPLITUDES 

In this section we describe a nondeterministic algorithm to prepare quantum 
states whose amplitudes are nonlinear functions of those of some input quantum 
state. 

Suppose we want to effect a quadratic transformation on the probability ampli- 
tudes Zj of \(f>). (We assume throughout that the initial state \<j>) can be efficiently 
prepared on a quantum computer, eg., using a combination of the methods of [7] 
and [8].) Since the amplitudes are unknown this is, in general, impossible and, 
unfortunately, this is the generic setting when integrating ODEs as every point on 
the solution trajectory can be regarded as an initial condition for the system. But 
suppose we have two copies of \<f>). In this case the probability amplitudes of the 
tensor product are given by 



(3) \m) = 2 E z Mjk), 

j,k=0 
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where, for convenience, we set zq = 1 from now on. Evidently every monomial 
z - z 1 ^ , Ij , Ik < 1, of degree less than 2 appears (more than once) in this expansion. 
Suppose we want to iterate the transformation 

(4) z^F(z), 
where 

/ 2 (z) 

(5) F(z) = . 
and / a , a = 1, 2, . . . ,n, are quadratic polynomials 

n 

(6) /a(z) = 51 a tl Z ^U 

kd=0 

with aj^ = ajj^ and /o(z) = 1- Thus we aim to prepare the quantum state 

1 " 

(7) i0')-- 7! E^( z )i a )- 

* a=0 

where, to simplify matters, we've assumed that the transformation is measure pre- 
serving, i.e., 

n n 

(8) l = ^|z,| 2 = ^/:(z)/ Q (z). 

j=l a=l 

The measure-preservation assumption plays an important role in our constructions, 
however, if it is relaxed our algorithm still proceeds unchanged: only the success 
probability is modified. 

To ensure that our method is efficient we need to make several extra assumptions 
beyond measure preservation. The first assumption is that \aff\ = O(l), k,l,a = 
1, 2, . . . , n. The second assumption is that the map F is sparse, which means that 

(9) |{(M)|aj£V0}|<*/2, a = l,2,...,n, 
and 

(10) \{<*\a$ ^0}| <s/2, k,l = 1,2,..., n, 

where s = 0(1). Note that the assumption of sparsity means that each / Q (z) can 
only involve at most s/2 monomials and that each variable appears in at most s/2 
polynomials / Q (z). The final assumption is that the Lipschitz constant for our 
system, i.e., that number A such that 

(11) ||F(x-y)|| <A||x-y|| 

in the ball ||x|| 2 < 1 and ||y|| 2 < 1, is 0(1). While these assumptions are rather 
restrictive, as we discuss, there are still many important systems which satisfy them. 
We also later describe how to relax these assumptions. 

It turns out that implementing the desired transformation will, in general, require 
that we make some destructive nonunitary transformation of the system's state. To 
understand what is required we now set up the operator 

n 

(12) A = H o$\aQ)(kl\. 

a,fc,Z=0 
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We now adjoin a qubit "pointer" P and use A to set up a hamiltonian (this is 
essentially the von Neumann measurement prescription [16]): 

(13) H = -iA®\l) P (0\+iA*®\0)p(l\. 
We now initialise our system in the state 

(14) mmp, 

and evolve according to H for a time t = e. The time we can evolve for depends 
crucially on the sparsity of H and the desired error [3] ; roughly speaking, when the 
sparsity of A is some constant s we can efficiently simulate (in terms of log(n)), up 
to some prespecified precision, the evolution for any constant time t. 
After the evolution the system ends up in the state 



I*) = e ieH \<P)\<P)\Q) = ^P-I^>l°) 



'-"mm = T, . 

(15) J=0 3 

= mm + e ^Mii)-.... 

Now noting that 

1 ™ 1 

(16) A\<f>)\4>) = - J2 a$z k z l \a)\0) = -=\<f>')\0) 

we now measure on the ancilla qubit and postselect on "1"; we succeed with 
probability w |e 2 (thanks to the measure preservation property; if our polynomial 
map is not measure preserving then this success probability will be proportionally 
lower) and our posterior state is 



(17) ^(i®i®<i|)|tt) = _L J2 al?z k z l \a)\0) = \cj>')\0). 

V 2 a,kd=0 

If we fail then we end up with some "poisoned" state (this occurs with proba- 
bility w 1 — |e 2 ), which we discard. In order to ensure, with high probability, that 
we end up with at least one copy of \<j>'} we need to repeat this process on roughly 
16/e 2 fresh pairs It is an interesting question whether one can design H so 

that the poisoned state can be recovered and used again. Such a possibility would 
allow one to substantially reduce the resource requirements of our algorithm. 

(Obviously, because the expansion (15) is truncated to first order, we don't 
exactly produce \<f>'), but rather some approximation to | </>'). To correct this we 
actually use the method described in [8] to implement the transformation 



(is) i^}|o> i ► ^i-e*m\4>)\4>m+^H\mm, 

where |e| < l/||i/||. Since H is sparse we have, by Gersgorin's theorem [10], that 
\\H\\ < cs, where c is a constant. This allows us to assume that when we succeed 
we will obtain precisely a copy of \<fi'}, modulo only imperfections in the simulation 
of e teH . We discuss these errors in the appendix.) 

If we want to iterate the polynomial map z F(z) a constant (in n) number m 
times to produce the state \^ m >) we need to start with (-jf ) m initial states. Thus 
the total spatial resources required by this algorithm scale as (If ) m log(n) because 
the cost of storing n variables via encoding in \<j>) scales linearly with log(n). (See 
Proposition 1 for a precise statement of the spatial resource requirements of our 
algorithm.) 

The simulation of the evolution e ltH on a quantum computer cannot be done 
perfectly; we quantify, in proposition 2, the errors that accumulate throughout 
the running of the polynomial iteration algorithm: if the evolution e' ltH can be 
simulated up to error S, then after m steps the final state will have accumulated 
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an error no worse that S(3 , y) m+ . Thus, choosing the simulation error to satisfy 
S < (3j)~ m will ensure that the mth iterate is exponentially close the desired state. 
The costs [3] of simulation to this level of precision imply that the total running 
time T of our algorithm scales as T <~ mpoly(log(n) log*(n))s 2 9 Kv/ ™, where n is 
some 0(1) constant and 

(19) log*(n) = min{r | log^ r) (n) < 2} 

with log^ denoting the rth iterated logarithm, and we've assumed that the trans- 
formation (18) is carried out in parallel on each of the remaining pairs in each 
iteration. Thus the parallelised temporal scaling is subexponential in m and poly- 
nomial in log(n). 

To complete the description of our iteration algorithm we need to describe how 
to read out information about the solution. After m iterations the system will 
be, up to some prespecified error, in the quantum state \4>^ m ">) encoding the mth 
iterate of F in the probability amplitudes. To access information about the so- 
lution we need to make measurements of the system. In principle, any hermitian 
observable M — Y^j k=o Mj,k\j)(k\ may be measured to extract information. Via 
Hocffding's inequality we learn that we can estimate, using repeated measurements, 
the quantity 

n 

(20) (M) = ZjMj,kZ k 

j,k=0 

to within any desired additive error. In practice the observable M is measured 
using (a discretisation of) von Neumann's measurement prescription [4, 16], so the 
evolution e ltM must be efficiently simulable on a quantum computer. Many natural 
such operators fall into this class, including - via a quantum fourier transform - 
the operator whose measurement statistics provide information on the sums 

(21) S fe = -=]>> e ^. 

As noted in [8] , and evident from our construction here, by measurement of multiple 
copies of \(t>^) it is also possible to extract information about polynomial functions 
of the solution. 

The flexibility to measure efficiently implcmcntable hermitian operators to ex- 
tract information about the solution provides the key to the exponential separation 
between our method and the best classical method. Indeed, if we only wanted to 
learn one element Zj , for some j, of the mth iterate then there is actually an efficient 
classical algorithm with the same resource scaling. (This is similar to the situation 
with linear equations [6, 2] where if we only want to learn about one element of the 
solution vector of a well-conditioned sparse set of linear equations we can do this 
efficiently classically. Indeed, even if the system is badly conditioned, we can still 
learn about parts of the solution in the well-conditioned subspace.) 

4. Solving nonlinear differential equations 

In this section we show how to use the method we've just described to integrate 
a sparse set of simultaneous nonlinear differential equations for any constant time. 

Suppose we want to integrate the system 1 (1) with the initial condition z(0) = b, 
where fj are sparse polynomials. The simplest approach is to use Euler's method 



We again assume that our system is measure preserving which means that X)j=0 z j /j( z W) + 
Zjf*(z(t)) = and that the Lipschitz constant is O(l). 
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[17]: we pick some small step size h and iterate the map 

(22) zj i ^ zj + hz'j = zj + hfj(z). 

While Euler's method is pretty terrible in practice, especially for stiff ODEs, it 
does provide a basic proof of principle; more sophisticated methods such as 4th 
order Rungc-Kutta and predictor-corrector methods are essentially only polyno- 
mially more efficient. As our algorithm is precisely an implementation of Euler's 
method on the probability amplitudes, we can appeal to the standard theory (see, 
eg., [11]) concerning its correctness and complexity; thus our algorithm suffers from 
all the standard drawbacks of Euler's method, including a lst-order decrease in er- 
ror in terms of the step size h. Nonetheless, it will provide us with an algorithm 
scaling polynomially with log(n), where n is the number of variables. However, it 
should be noted that our method scales exponentially with the inverse step size. 
Thus, without modification, our algorithm is really only suited to well-conditioned 
systems. 

As we've indicated, the idea behind our approach is very simple: we integrate 
the system using Euler's method. Thus, given [</>(£)) we aim to prepare 

(23) m+h)) = \m+h\4>'(t))+o(h 2 ), 

where now 

1 n 1 n 

(24) |0'(t)> = -^E/aW*))li> = -^ E a$z k {t) Zl {t)\a). 

VZ a=0 VZ Q ,fc,Z=0 

To implement this transformation we suppose we have two copies of \<fi(t)) and apply 
the method of the previous section to implement the polynomial transformation 
z a !— » z a + hf a (z(t)). Note that this transformation is only measure preserving to 
0(h); the success probability will be diminished by a factor of 0(h 2 ), which can be 
made negligible by reducing h in the standard way. 

So, to integrate the system (1) forward in time to t = 0(1) we begin by dis- 
cretising time into m steps. (Thus our step size is h — t/m.) We then prepare 
(If)" 1 copies of |^(0)), where e is as in §3, and apply the method of §3 to produce 
approximately (if)" 1 copies of \(f>(t/m)) in expected time poly(log(n)). We then 
iterate until we produce at least one copy of \4>(t)) in expected time poly(m, log(n)) 
with probability greater than 1/3. The resources required by this approach scale 
polynomially with log(n) and exponentially with t and l/h. 



5. Extensions and applications 

In this section we briefly describe several extensions and applications of our 
quantum Euler's method. 

What sort of systems will be tractable with our approach? We only sketch 
a couple of examples here, leaving the wider application of our approach to more 
detailed investigations. Because the sparsity and the measure-preserving properties 
play a key role in the resource scaling of our algorithm it is desirable to focus on 
those sparse systems preserving the "hamiltonian" Y^j=i \ z j\ 2 ■ O ne example of such 
a system is the Orszag-McLaughlin dynamical system [14, 15]: 

(25) ^ = Xj+ixj+2 +xj-iXj-2 - 2xj+iXj-i, j = l,...,n, 

with periodic boundary conditions x n+ \ = x\. The variables Xj are real and pre- 
serve Y^j=i x j- The dynamics generated by this system are extremely complicated. 

Another example of a system which can be studied using our algorithm is the (dis- 
crete) nonlinear Schrddinger equation on any finite graph G = (V, E) of bounded 
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degree: 



(26) 



.dZy 



2deg(v)z v - ^2 



v e V, 




w 



'V 



where k e N. 

Several extensions of our algorithm are possible. The first obvious extension is 
to systems whose nonlinearity is cubic or higher. This can be done in the natural 
way by consuming 3 (or more, for higher degrees of nonlinearity) copies of \4>(t)} at 
each step to implement the desired nonlinear transformation. A second extension 
is to certain densely defined systems, i.e., those for which the A operator is dense: 
because, following [8], we can actually implement any efficiently computable func- 
tion g of A in the step, we can access some dense operators g{A). A third extension 
allows the efficient computation of the equal-time statistics of deterministic dynam- 
ical systems. Here the idea is to apply our method not to \cf>} = z a (t)\a), but 
rather to one half of an initial state which is an entangled pair: 



where (z|z') = 5(z — z') and d//(z) is an efficiently implementable probability mea- 
sure, eg., uniform or gaussian [7]. Applying our algorithm to this initial state allows 
us to efficiently sample equal-time statistics via measurements on P. This should be 
contrasted with the classical Hopf functional approach [9] to solving this problem 
which introduces a Fokker-Planck type partial differential equation to study distri- 
butions of solution trajectories. (The application of the Hopf functional approach 
to the Orszag-McLaughlin system is considered in [12].) 

Finally, it is not implausible that there is a trade-off between time and space; per- 
haps there is a quantum algorithm which integrates a constant number of variables 
which scales polynomially with log(t) and — log(/i)? 



We have presented a quantum algorithm to iterate large systems of sparse polyno- 
mial maps. We've also described an implementation of Euler's method to integrate 
a system of ODEs. As long as the system is sparse the resources required by the 
method are polynomial in log(n), where n is the number of variables. However, the 
resources consumed by the method scale exponentially with the inverse step size 
and the integration time, as well as with the degree of the nonlinearity. 
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Appendix A. Proofs of the claims 

In this appendix we present the technical proofs of the claims made in the text. 

A.l. Description of the algorithm. In this subsection we provide the formal 
specification of our iteration algorithm. 



Algorithm 1 

Set N — (^)~ m , according to Proposition 1; 

Initialise the system in the state |$) := (|</>}|</>)|0)) l8,JV / 2 , where \cf>) = ^|0) + 

Set H := -iA® |1)(0| +iA* ® |0)(1|; 
for i = m to 1 in steps of —1 do 
5 = 0; 

for j — 1 to N in steps of 1 do 

Evolve the jth pair |0 2 j_i)|02j)|O) according to (18), measure the ancilla, 

and postselect on "1" . 

:= -iV2(I ® I ® (l|)tf I^-OI^IO); 

if success (i.e. |<^) = |#y-i>l°» then S:=S + 1; 
end for 
7V:=2LS/2J; 

if S < 2 1 ~ 1 then Exception: algorithm failed; 

Set |$) :=(|0')|0')|O)) 8JV / 2 ; 
end for 



A. 2. Bounding the expected running time. We now prove the following 

Proposition 1. Let m > 6. If the number N of initial states for Algorithm 1 
with p = e 2 /2 is (^) _m then it succeeds in producing at least one state |(/>( m ') with 
probability at least 1/3. 

Proof. Let S be the random variable which counts the number of successes out of 
the N/2 trials in each round and let Nj denote the number of successfully pro- 
duced states in round j. Hocffding's inequality provides a bound on the cumulative 
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distribution function of a binomial random variable S with parameters (N,p): 

2(Np-fc) 2 

F(k; N,p) = P(S <k)< e * 

for < k < E [S] = ^N. We declare failure if S < k, so that, defining k = XN, we 
have 

_ 2(Wp-AW) 2 

P(failure) < e « 

= e -2N(p-X) 2 ^ 

Therefore 

P(success) > 1 - e - 2JV (f" A ) 2 . 

We want the TO-step algorithm to produce at least one state \<fr( m >) encoding the 
mth iterate, i.e. we want N m > 1 with probability at least 1/3. That is, we require 

P(m successes) = P(Successi)P(Success2) • • • P(Success m ) 

> (1 - e- 2Nl{p - X)2 )(l - e - 2Ar 2(P- A ) 2 ) ... (1 - e -2N m (p-\) 2 ^ 

> 1/3. 

If the algorithm succeeds, after one step we produce at least | = states (A; 
successes) with probability at least (1 — e ~ 2Nlt - p ~ x ) 2 ), and so after j steps we have 
that 

< A x,-i 



(28) Nj >^-j N 

with probability at least (1 - e" 2iv i(P- A ) 2 )(l _ e -2N 2 ( P -x)^ ...(!_ e -2N 3 ( P -xfy 
To ensure that the final success probability is greater than 1/3 we demand that 



P(Successj-) > 1 - e - 2 ^(p- A ) 2 > l - 1 



m 
so that 

P(m successes) > ^1 - — ^ = e" 1 + 0(l/m) > 1/3, 
for to > 6. Therefore, for all j we need that 

l _ e -2N ] ( P -Xf >!_! 

TO 

that is, 

(29) - > e -2^(P-A) 2 . 

TO 

The RHS of (29) is decreasing in Nj and so is maximum when j — to. From (28) 
we have 

N m > (-) N 

substitution into (29) gives 

1 > e - 2 ur~ lAr (p- A ) 2 , 

TO 

that is, we require that 

/ \ \ m ~ 1 
logTO<2 - N(p-\) 2 . 
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If we choose A = p/2 and N = (f )" m = (|) _m then from (28) we have 



4- 

V,„ I - N 

A\ m - 1 /A 



2/ V 4 

ry rn ^ 

" A ' 

and the result follows. □ 

A. 3. Bounding the accumulated error. We now bound the accumulated error 
that builds up throughout the running of the quantum iteration algorithm. 

Proposition 2. The error S m — \\\4>^) — |^ m ')|| that accumulates after m iter- 
ations of the quantum iteration algorithm 1 is bounded by 

,/(3 7 r +1 -i 

m - 3 V 37-1 

where r\ is the error in simulating e leH , |^M m )) is the actual state produced by the 
algorithm, and 7 is an 0(1) constant which depends on the sparsity s. 

Proof. Our proof works by analysing the errors that accumulate in pairs 
of states during the jth round. Suppose we have an error in the starting pair, so 
instead of |0o) we have \(f)o) + |A<^o) = IV'o), say, where the initial error |A0 O ) has 
magnitude So, i.e. So = |||A0o)||- Suppose also that our simulation of U — e leH is 
imperfect, i.e. V is the operator that is actually applied and \\U — V||oo < V- 

We initialise our system into the state IV'o^V'ollO) and evolve according to V. 
This particular pair will then be in the state 

^o)|^o>|0) 

To measure this state on the ancilla qubit and postselect on "1" we apply the 
measurement operator Pi = I® I® |1)(1|- The (subnormalised) posterior state of 
the pair is then 

|^> = P 1 V(\<h) + |A^,))(|0o> + |A0 O ))|O). 

Recall that if there were no errors in the starting state and the algorithm was 
perfect we would have 

V2- 



} = — Pitf|0o>|0o>|O>. 



Let the error between the subnormalised posterior states |?/>i) and 1^) after one 
step have magnitude Si — \\\tl>i) — \4>i)\\- We bound this as follows 

Si < \\I\V(\A<M\<h) + \<h)\&<fo) + |A^o)|A0o)) + Pi(V - t/)|<MI<MH 
(30) < |||A0 o )|0o) + |^o)|A0 o ) + |A^ o )|A0o)!l + \\(V- U)^ 

< 5 + S Q + Si + 77 < 3<5 + r}. 
We bound the error Si between the normalised posterior states by 

= — —(3S + rf) = li'ASo + rf), 
where 7 is, by assumption, a constant of order s. 
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Repeating this argument allows us to set up the recurrence 

(32) Sj < 7(3^-i j = l,2,...,m. 

Solving the recurrence, and assuming that the initial states are constructed perfectly 
(i.e., S = 0) gives us 

(33) , m <_^____ij. 

□ 
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